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We address the crystallization of monodisperse hard spheres in terms of the properties of finite- 
size crystalline clusters. By means of large scale event-driven Molecular Dynamics simulations, we 
study systems at different packing fractions (p ranging from weakly supersaturated state points to 
glassy ones, covering different nucleation regimes. We find that such regimes also result in different 
properties of the crystalline clusters: compact clusters are formed in the classical-nucleation-theory 
regime (</> < 0.54), while a crossover to fractal, ramified clusters is encountered upon increasing 
packing fraction {(j) > 0.56), where nucleation is more spinodal-like. We draw an analogy between 
macroscopic crystallization of our clusters and percolation of attractive systems to provide ideas on 
how the packing fraction influences the final structure of the macroscopic crystals. In our previous 
work {Phys. Rev. Lett., 106, 2i570i, 20il), we have demonstrated how crystallization from a glass 
(at > 0.58) happens via a gradual (many-step) mechanism: in this paper we show how the 
mechanism of gradual growth seems to hold also in super-saturated systems just above freezing 
showing that static properties of clusters are not much affected by dynamics. 



I. INTRODUCTION 

Following the pioneering computer simulations of 
Alder and Wainwright [J , more than 50 years ago, assem- 
blies of hard spheres in thermal motion have become the 
subject of intense research aimed at unveiling the funda- 
mental physics of both thermodynamic and kinetic phase 
transitions. This numerical work has been complemented 
by experimental studies of suspensions of "hard-sphere" 
colloidal particles (see Ref. [2-5]). 

Despite the simplicity of the hard-core interactions, 
the behaviour of this system is far from trivial and sev- 
eral important questions still remain unanswered. Con- 
cerning thermodynamics, the location of the fluid-to- 
crystal transition, occurring in the packing fraction (f) 
range 0.492 < (p < 0.543 has been quite well charac- 
terised [Hlg-ll]. However, open issues remain concerning 
the effect of polydispersity |9l-[ll| and, most importantly, 
the reported discrepancy between calculated nucleation 
rates and those measured experimentally [12!]. These is- 
sues have recently caused a revival of studies on hard- 
sphere nucleation [isl-llq. In particular, recent works 
aim to elucidate the mechanism of crystal formation in 
terms of preferred precursor structures jlTHig ]. propos- 
ing a two-step scenario [20, [21|, and also aim to assess 
the competition between different solid polymorphs [23 ]. 

Not only does the fluid-to-crystal transition remain 
debated, but also the presence of an ideal glass tran- 
sition at packing fractions below random close packing 
((j^rcp ~ 0.64) is still under discussion. The initial obser- 
vations by Pusey and van Megen [2] of a glass transition 
taking place for ^~ 0.58 were later challenged by works 
in microgravity |23l . |24| and by numerical simulations by 
Torquato and coworkers [25|. Moreover, a recent work 
by Brambilla et al. [2g caused considerable debate in 



the glass community [27H29J . By considerably extending 
the observation time window, these authors suggested 
that the glass transition was always obviated by acti- 
vated processes whenever (jj < (jjrcp- However, delicate 
issues remain, concerning for example the correct way of 
estimating the "true" colloidal packing fraction in these 
experiments [30| . 

With this as background, we have recently started an 
extensive numerical investigation of hard-sphere systems, 
via event-driven Molecular Dynamics simulations. Ini- 
tially, we assessed the role of polydispersity on dynamics 
and crystallization [13|, |31|, finding that, while changing 
substantially the nucleation behaviour, a small degree 
of polydispersity s < 6% does not strongly affect the 
dynamics, leaving unchanged the tendency to form an 
ideal glass around 0.58 |31j. In addition, we identified 
two regimes of nucleation [13l . |33] : standard nucleation 
and growth at concentrations in and slightly above the 
coexistence region and, beyond (j) '^ 0.56, "spinodal nu- 
cleation" , where the free energy barrier to nucleation ap- 
pears to be negligible. In the latter regime, at very high 
volume fractions, we have investigated the subtle inter- 
play between slow dynamics and crystallization, eluci- 
dating a novel mechanism by which hard-sphere glasses 
crystallize without particle diffusion (even at the single 
particle level) on scales comparable to or beyond their 
radius [33| . Diffusionless crystallisation has also been re- 
cently reported for a 3 dimensional lattice system [3J]. 

In this work, we investigate an aspect of crystallisa- 
tion of monodisperse hard spheres that, to our knowl- 
edge, has not been studied so far: the statistics proper- 
ties of crystalline clusters as crystallisation proceeds. By 
monitoring a large system, we identify the solid-like clus- 
ters and analyse their size distribution and shape as a 
function of packing fraction and of time. These clusters 
undergo a percolation transition of the kind described 



in simulation [35| and in theoretical work on molecular 
systems [3g|. We monitor the cluster size distribution 
and clusters radius of gyration as a function of time, to 
address whether non-trivial exponents can be extracted 
and related to the standard percolation scenario that is 
commonly encountered in attractive systems [37| . 

A comment is due on the definition of clusters of hard- 
sphere particles, where no attraction is present. Indeed, 
clusters are rigorously defined [3^ for attractive systems 
in terms of an energy scale which cohesively keep par- 
ticles together. Moreover, in the case of colloidal parti- 
cles interacting with competing short-ranged attraction 
(e.g. depletion) and long-ranged repulsion (induced by 
electrostatics), clusters are found to be stable minima of 
the underlying potential [39|-|42| . These equilibrium clus- 
ters have been observed in both experiments and sim- 
ulations |43l - t4a | and their role as a building blocks of 
non-equilibrium structures like Wigner glasses and gels 
is still under active investigation [40, |47h52i] . In the case 
of hard spheres, particles forming a cluster "glue" due to 
the thermodynamic driving force toward crystal nucle- 
ation |53l . |54| rather than due to attractive interactions 
between particles. Hence, whenever we find two nearest- 
neighbour solid-like particles (according to the definitions 
discussed below), we consider them as belonging to the 
same crystalline cluster. Of course, these clusters can be 
either transient or permanent, in analogy with attractive 
clusters at finite temperature. 

Within our study, we cannot monitor the equilibrium 
distribution of such clusters because of the irreversible 
process of crystallization once it is triggered. However, 
we can follow their growth while crystallization proceeds, 
so that with time clusters become larger and larger and 
consequently more and more permanent, until macro- 
scopic crystallization of the samples occur. In particular, 
we study crystallization at different packing fractions, in 
order to address the question of whether significant dif- 
ferences can be found between the different nucleation 
regimes. Our findings suggest an important crossover 
from compact to fractal clusters with increasing super- 
saturation, which reflects the change in the underlying 
crystallisation processes. These findings could have im- 
portant consequences for the final structures of the re- 
sulting macroscopic crystals. 

The manuscript is organised as follows: we first present 
the simulation and analysis method in Section [Hi next, 
we report our results on the clusters structure at different 
packing fractions (Section IIIIA[) . the percolation of the 
crystalline clusters (Section |IIIB|) . the cluster size distri- 
bution (Section IIII Cp and on the effect of the criterion 
to identify solid-like particle on properties of the growing 
crystalline clusters (Section Fill Pp . 



II. SIMULATION AND ANALYSIS METHODS 

We perform event-driven Molecular Dynamics simula- 
tions in the NVT ensemble with cubic periodic boundary 



conditions for a large system of N = 86400 monodisperse 
hard spheres [5a, [5g|. Mass, length, and time are mea- 
sured in units of particle mass to, particle diameter a 
and yjmcr'^ / kbT, where k,b is the Boltzmann constant 
and T the temperature and we set kbT — 1. The system 
is prepared at different packing fractions (f) — ^Na^ /V 
(with V the system's volume) beyond freezing ranging 
from (j) = 0.54 to = 0.61. At </) > 0.54, the metastable 
fluid-phase spontaneously crystallizes within the dura- 
tion of our simulations. With increasing 0, the relaxation 
of the system becomes slower and slower, until the glass 
transition around ~ 0.58 is encountered [3l|. Above 
this value, ageing effects are present and crystallization 
occurs from a glass [SSj . 

To study the formation and growth of crystalline clus- 
ters, we first identify the solid particles in the system. 
To this end, we use the rotationally invariant local bond 
order parameter dg defined as the scalar product between 
particle i's (jg complex vector and the one of each of its 
neighbour J, dG{i,j) (see Appendix Elfor further details). 
Once we have identified all solid particles Ng in the sys- 
tem, we monitor the fraction of crystalline particles, de- 
fined as X = -5^. To identify crystalline clusters, we run 
a cluster-algorithm: starting from a solid particle, we lo- 
cate its solid neighbours within a distance of 1.4cr and 
define them as belonging to the same cluster. After hav- 
ing iteratively applied the cluster-algorithm to all solid 
particles in the system, at every time-step we identify all 
crystalline clusters and compute their size (s). 

To improve the statistics of our results, we consider 10 
independent crystallisation trajectories for each studied 
packing fraction, each initiated from configurations with- 
out pre-existing crystal nuclei. In order to achieve this 
non-trivial condition, we compress a small system of 400 
particles to a high packing fraction, (j) = 0.64, monitor- 
ing that the fraction of crystalline particles is less than 
0.005 of the total number of particles. We then replicate 
this system periodically in space, checking that any sign 
of periodicity resulting from the replication procedure is 
lost after a very short time (as in Ref. [3J|). Next, be- 
fore starting Molecular Dynamics runs, we isotropically 
expand the configuration to the desired packing fraction. 

First of all, we perform an analysis of the structure of 
the growing crystalline clusters at different packing frac- 
tions. This is is done both qualitatively, monitoring the 
structure of the clusters that are visible in the snapshots 
taken from our simulations and quantitatively, studying 
the size-dependence of the cluster's radius of gyration: 
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where s is the number of particles in the cluster (i.e. 
the cluster size), fi indicates the position of particle i 
and RcM is the position of the centre of mass of the 
cluster. For fractal aggregates, the dependence of the 
radius of gyration on the cluster size follows a power 



law, i.e. Ra 
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, where df is the fractal dimen- 



sion of the aggregation process. For random percola- 
tion in three dimensions, the value of df is known to be 
about 2.5 [33] • For spherical clusters Rg '^ s^/"^, whereas 
for planar (loose) clusters Rg ~ s^^^ and for linear ones 
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Next, we study the percolation of crystalline clusters 
by estimating the first time when a system-spanning clus- 
ter appears in the simulation box, i.e., the percolation 
time Tp. In practice, we detect Tp as the time when the 
largest crystalline cluster becomes infinite according to 
our periodic boundary conditions criterion. After having 
identified, at each time-step, the largest crystalline clus- 
ter, we replicate the simulation box and check whether 
we get a single cluster in the replicated system: the time 
when this happens defines Tp (which is then averaged over 
independent runs). Once Tp is known, we estimate the 
value of the total crystallinity at percolation X{Tp) — Xp. 
Even though we are aware that an accurate location of 
the percolation threshold would require a finite-size scal- 
ing analysis, our aim is not to define with high precision 
the percolation threshold, but rather to find a link be- 
tween the percolation transition and macroscopic crys- 
tallization. 

We also calculate the cluster size distribution, i.e. the 
fraction of clusters with s solid-like particles. The clus- 
ter size distribution changes with time during the crys- 
tallization process. To reduce its numerical noise, we 
evaluate it for configurations of independent trajectories 
having approximately the same crystallinity, using nar- 
row X intervals of width AX = 0.01. In this way we can 
average our results over all trajectories within a given 
X-interval and over the 10 independent runs. According 
to random percolation theory, the cluster size distribu- 
tion in a three dimensional system decays exponentially 
fast for large s below the percolation threshold, whereas 
it follows a power-law behaviour at percolation, i.e. s~'^ 
(where t ~ 2.2 [13). 

Finally, to assess the role of the particular choice of pa- 
rameters used to define solid-like clusters, we analyse the 
effect of the "solid-like" criterion on the clusters radius 
of gyration, considering different choices of the number 
of solid-like connections £,c needed to identify a solid-like 
particle. 



III. RESULTS 



We present our results as follows: in Section fill Al we 
perform both a qualitative and a quantitative analysis 
of the clusters' structure for packing fractions ranging 
from <j) — 0.54 to ^ = 0.61; in Section Fill Bl we study 
percolation of the crystalline clusters; in Section llllCl we 
estimate the cluster size distribution. In the last Sec- 
tion, IIIIDl we analyse the effect of the criterion used to 
identify solid-like particle on properties of the growing 
crystalline clusters. 



A. Structure of the growing crystalline clusters 

First, we study the characteristic structure of grow- 
ing crystallites at different packing fractions (p. In or- 
der to visualise how the structure of crystalline clusters 
changes with increasing packing fraction, we follow the 
largest cluster during the simulation run, as represented 
in Figure[T] When (f> = 0.54 (left-most panel in Figure[T]), 
we observe a rather homogeneous and compact cluster, 
whereas more ramified structures emerge at higher (j) (e.g. 
(j> — 0.56,0.58,0.61). An intermediate case is = 0.55 
(second panel from the left in Figure [IJ clearly showing 
small compact clusters connected to each other via thin 
branches. 

As we have already discussed elsewhere [ij, Isj] , there 
is a strong correlation between the crystallization mech- 
anism and the packing fraction of the system. For (j) = 
0.54, crystallization appears to follow the classical nucle- 
ation theory (CNT) picture where a critical cluster needs 
to form in order for crystallization to proceed. Once a 
crystalline cluster overcomes a critical size, it grows in an 
almost compact structure and the entire system crystal- 
lizes ( "classical nucleation and growth" regime) . Having 
used a large system and being at a relatively high su- 
persaturation (with respect to freezing), we detect more 
than one critical cluster in the system (of the order of 2- 
3): all these clusters are compact, can form anywhere in 
the system, grow independently and, eventually, merge 
together. When cj) > 0.56 ( "spinodal-like crystallization" 
regime), a random crystal growth habit is reached, where 
crystallization takes place everywhere in the system. The 
free-energy barrier of nucleation gets closer to zero the 
higher the packing fraction, and no critical cluster size 
needs to be exceeded in order for crystallization to pro- 
ceed. Crystallization proceeds via branching of ramified 
clusters as X grows (right-most snapshot in Figure[T]) and 
this mechanism is reminiscent of a percolation transition 
associated with the crystallization process: moreover, at 
(j) > 0.58, crystallization happens without particles dif- 
fusing more than a diameter [SJl ■ As said above, the case 
at (j) = 0.55, where branching connections of more com- 
pact (and finite) clusters are found, is intermediate be- 
tween (p — 0.54 and higher packing fractions {4> > 0.56). 

To quantify the structure of the crystalline clusters, 
we measure their radius of gyration (Rg) as a function 
of their sizes at all packing fractions and extract from 
it an effective fractal dimension df, via Rg ~ s^^'^f (see 
Figure m. To improve the statistics, data are averaged 
over all equal-sized clusters, over 10 trajectories and over 
time. 

In Figure [2] we observe that, at all packing fractions, 
the radius of gyration for small clusters follows Rg ^ s^l"^ . 
When clusters become larger (s > 50), their fractal di- 
mension starts to be larger than 2 and their internal 
structure changes depending on 0. Two distinct be- 
haviours emerge: one pertaining Xo p — 0.54 and an- 
other to > 0.56 (as shown by the superposition of the 
data in Figure [5]). When = 0.54, crystallization fol- 
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FIG. 1: Snapshots of typical largest crystalline clusters as a function of at fixed cluster size (s = 5000). This size is achieved 
when {(l>,X) are respectively (0.54, 0.08), (0.55, 0.15), (0.56, 0.15), (0.58, 0.14), (0.61,0.12). Periodic boundary conditions are 
taken into account and clusters are centred in the simulation box. 
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FIG. 2: Radius of gyration {Rg) versus clusters-size (s) at 
various (j> (s is averaged over different runs and over time). 
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As a guide to the eye, we plot Rg 

(continuous black curve), 2.5 (dashed magenta curve) and 3.0 

(dashed blue curve). 



cess any more (there is no free-energy barrier to over- 
come [13| ) and all Rg data scale onto the same master 
curve: after a loose/open growth of the small clusters, 
as soon as s > 50, clusters grow following the scaling 
predicted for random percolation in 3D with fractal di- 
mension df '^ 2.5. This is consistent with the snapshots 
presented in Figure [1] where already at (f> — 0.56 clus- 
ters start to be ramified and crystallization occurs with 
the growth of clusters appearing almost everywhere in 
the system. However, as we will discuss in Section UlI CI 
clusters do not grow completely at random; instead they 
tend to form in the vicinity of already existing ones [33| . 
This is even more evident at packing fractions beyond 
(j) — 0.58, due to particles' lack of diffusion on the scale 
of their radius. 

Finally we discuss the case (f> = 0.55, where we ob- 
serve signatures of an intermediate behaviour. Given 
that the nucleation free-energy barrier is quite low, clus- 
ters start growing everywhere in the system, create a 
branch-like structure when they grow larger (as shown in 
the snapshot of Figure [T]) , and Rg shows an intermedi- 
ate behaviour between the compact scaling of the CNT 
regime and the fractal one pertaining to the spinodal-like 
regime (as shown in Figure [5]). 



lows the standard CNT picture, i.e. the system has to 
overcome a free-energy barrier in order to transform into 
a crystal. Hence, nucleation is a rare event and clusters 
appear and disappear in a stochastic way until the crit- 
ical size is reached. The critical cluster size at this (j) i s 
estimated with CNT to be of the order of Sc '^ 50 \13^ . 
As it can be observed from Figure [H beyond this criti- 
cal value clusters show a marked change in their fractal 
dimension, becoming more and more compact (value of 
df increasing). When the cluster size exceeds « 300, the 
data clearly follow the exponent d/ = 3, characteristic of 
a compact cluster (such as the corresponding cluster in 
Figure[T]) [S^. Hence, we conclude that, at this c^ (and in 
the CNT-regime) , the structure of the clusters is mainly 
determined by their size. 

When (f) > 0.56, nucleation is not an activated pro- 



B. Percolation analysis 

We now perform an analysis of the percolation of the 
largest crystalline cluster in the system for all studied 
packing fractions, identifying both the percolation time 
Tp and the value of the crystallinity at percolation Xp. 

In Figure |3] we report the dependence of Xp on and 
observe that, for increasing values of (/>, the crystallinity 
at the percolation transition has a monotonically decreas- 
ing behaviour, that is steeper for < 0.56 than for 
(/) > 0.56. To interpret this result we recall our previous 
findings that clusters are more compact with decreasing 
(/) for < 0.56. Hence, in order to percolate, the largest 
cluster needs to reach a larger size than the one needed 
at higher (j>. On the other hand, when (j> > 0.56 df re- 
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FIG. 3: Fraction of crystalline particles at the percolation 
threshold Xp as a function of (j). Inset: percolation time Tp as 
a function of A. 



mains constant so that Xp does not change much with (j). 
Due to the stochastic nature of nucleation at (j> = 0.54, 
that results in a small number of growing clusters that 
is different for each independent run, the Xp error bar is 
larger at the lowest packing fractions. 

For (f> > 0.56 crystalline clusters are found to perco- 
late when the total crystallinity of the system is Xp sa 
0.10 — 0.15. The fact that the overall crystallization needs 
to be at least 0.1 in order to observe percolation of the 
crystalline clusters (made of purely repulsive spheres) has 
an analog in studies of percolation of attractive clusters 
(made of particles interacting via short-range attraction 
and screened electrostatic repulsion) |45j|. In fact, it has 
been shown that in such systems, when the colloidal 
packing fraction 0c is about 0.10, attractive clusters 
(whose size distribution obeys a power law with an expo- 
nent typical of random percolation) percolate throughout 
the system. 

So far, we have presented only static observations of 
the structure of the growing clusters. In the inset of 
Figure [3] we plot kinetic results of the clusters growth, 
representing the percolation time Tp as a function of 0. 



Contrary to the monotonic decrease of Xp, Tp displays a 
clear minimum at intermediate 0. The fastest percola- 
tion of the largest cluster is at </> = 0.56 when the largest 
cluster is already branched but particles can still easily 
diffuse. Whereas at low cj) the percolation time is longer 
due to the compact structure of the growing cluster (re- 
quiring a larger Xp), at large percolation slows down 
due to the slowing of particles dynamics. This behaviour 
is remarkably similar to that observed for the "nucleation 
time" Tn{(f>), defined as time at which X — 0.2, that cor- 
responds to the time when the system is strongly commit- 
ted to crystallize [1J|: at low 0, nucleation is slow (large 
T„) due to a high nucleation free-energy barrier, while at 
large (j) nucleation slows down again due to the slowing 
particle dynamics; the fastest nucleation (smallest r„) 
occurs around (j) = 0.56, where the free-energy barrier is 



very low and diffusion still possible. Therefore, the per- 
colation of the growing clusters is correlated to the com- 
mitment of the system to fully crystallize. We stress that 
at large 0, despite the slowing down of particles dynam- 
ics, percolation is achieved with lower crystallinity than 
at lower 0. Hence, the kinetic slowing down is probably 
responsible for the continuous decrease of Xp, even when 
the shape of the clusters (quantified by df) stops chang- 
ing. Indeed, we find that Xp{(j) = 0.61) is almost a factor 
of 3 smaller than Xp{(j) — 0.54), whereas Tp{(j> = 0.61) 
is around 3 times larger than Tp(0 = 0.54), due to the 
slowing down of the dynamics beyond the glass transition 
(occurring at ~ 0.58 [i,[3H). 

The behaviour of Xp with (f) can also be rationalised by 
looking at the total number of clusters Nc- In the inset 
of Figure m we report Nc as a function of X for different 
packing fractions. 
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FIG. 4: Xp (black curve) and Nc{X = 0.10) (red curve) as a 
function of 0. Inset: Nc as a function of X (averaged within 
a AX = 0.01 interval and over independent trajectories) at 
different (j) (from bottom to top curve: =0.54, 0.55, 056, 
0.58, 0.61.) 



It is evident that, for a given X , the number of clus- 
ters increases with increasing 0, due to the reduced free- 
energy barrier, and it saturates when > 0.58. In partic- 
ular, for (p = 0.54, the maximum number of clusters oc- 
curs at very small values of X, where one or more critical 
nuclei have already formed and grow by single particles 
or small clusters attaching to them (that explains the 
slow decrease of Nc with increasing X). On the other 
hand, with increasing 0, Nc keeps increasing for small 
values of X (since many small clusters can easily form) 
until it starts decreasing from X ^ 0.03 — 0.05 (where 
clusters start to branch). We also observe a clear depen- 
dence on (j> below the glass transition (0 = 0.58). This 
can be better seen by looking at Nc for a fixed value of 
X, e.g. X — 0.10, which is reported in Figure 0] together 
with Xp as a function of (f). Comparing the behaviour of 
Nc to Xp, we can rationalise the decrease of the latter as 
being related to the increase of the total number of clus- 
ters up to 4> ^ 0.58. For larger cj), Nc remains constant 
and the decrease of Xp could be partially ascribed to the 
increase of packing fraction, so that even if the number 



of the clusters is the same, the number of chisters per 
unit volume is higher and percolation happens at lower 
X. Moreover, minor differences within the statistical er- 
ror of our analysis of the clusters shape and the cluster 
size distribution (see next Section) could also contribute 
to the decrease of Xp from (j) — 0.58 to = 0.61. 



C. Cluster size distribution 

In this section we discuss the behaviour of the cluster 
size distribution. n{s) represents the number of clusters 
of size s and /(s) the fraction of such clusters, i.e. n{s) 
divided by the total number of clusters N,,. n{s) varies 
as crystallization proceeds. To improve the statistics, we 
average over 10 independent trajectories and within X 
intervals of width 0.01. For all packing fractions we cal- 
culate n{s) and /(s) for crystallinity up to X = 0.15; at 
this value, crystalline clusters have not percolated yet at 
(j) = 0.54 and (j) = 0.55 (Figure |3]) and have just perco- 
lated (on average) for > 0.56. However, in the analysis 
that follows we have excluded percolating clusters. 

In Figure [5] we report n{s) for different X-intervals 
both for (j) = 0.54 (top) and = 0.61 (bottom). For 
both packing fractions, n{s) shows a rapid decay for 
X < 0.01 signaling the growth of the first few crys- 
talline clusters. Once X > 0.03 (when (j) = 0.54) and 
X > 0.06 (when (j) = 0.61), all curves collapse onto each 
other, showing that the clusters distribution has reached 
a stationary profile during crystal growth. This hap- 
pens for all studied (j) and justifies further averaging of 
the n(s) curves over all configurations within the range 
0.06 < X < 0.14 [H]. We also notice that large clus- 
ters, e.g. s > 100, are rare in the case oi (j) = 0.54 
while their number becomes very large for = 0.61, as 
indicated respectively by scarce/abundant population of 
large-s points in Figure. [5] 

We now focus on the fraction of clusters of size s, /(s). 
In order to evaluate /(s) for increasingly large s (where 
clusters of each size become more and more scarce), we 
apply the following criterion: (i) we arbitrarily subdivide 
the range of s- values into intervals Asi = {si^max — Si^min) 
containing at least one cluster each and estimate n*, i.e. 
the total number of clusters within each interval Asf, (ii) 
we assign the value of rii — n^/Asi to every s within 
Asi. The total number of clusters can be computed as 
Nc = X)i "-i^^i s-nd the fraction of clusters of size s, 
/(s) = Jii/Nc- We now represent /(s) as a function of s, 
the closest integer to the central value of As^. We leave 
to the Appendix IB] the definition of all Asi — {si^max — 
Si,min) intervals. In Figure H] we plot /(s) for several 
packing fractions ranging from 0.54 to 0.61. Data are 
averaged for 0.06 < X < 0.14 (within this range we have 
found the profiles to be stationary and independent of 

X H). 

First of all, we notice the marked difference that is 
found between the data corresponding to cf) — 0.54 and 
all other data. This demonstrates that the stationary 
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FIG. 5: Number of clusters riiyS] at = 0.54 (top) and c^) — 
0.61 (bottom). Each data set is averaged over 10 independent 
runs in intervals ranging from 0.0 < X < 0.01 to 0.14 < X < 
0.15 (AX = 0.01). n{s) can take values smaller than 1, being 
computed as the number of clusters for each value of s divided 
by the number of runs and by the number of configurations 
analysed for a given X window. Symbols are the same for 
both panels. Only non-percolating clusters are considered in 
the analysis. 



cluster size distribution is, within our statistical uncer- 
tainty, identical for any larger than 0.55. This be- 
haviour is similar to what we observed for Rg, pointing 
to the universality of (static) properties of the clusters 
once spinodal-like nucleation has occurred. 

From the calculated /(s) we can search for the emer- 
gence of a power-law behaviour approaching crystalliza- 
tion, in analogy with standard attractive systems ap- 
proaching a percolation transition. We find that when 
(f) = 0.54 a crossover between two regimes takes place at 
s* = 50, which corresponds to the estimate of the crit- 
ical cluster size according to CNT [1J|. For s < s* the 
cluster distribution is not too different from that of a 
system approaching random percolation (/(s) '^ s^^'^); 
whereas larger clusters grow in a compact way, result- 
ing in a much lower power-law exponent (/(s) '^ s"^'^). 
When (p > 0.55, f{s) seems to follow a power law with 




FIG. 6: Fraction of clusters of size s, /(s), averaged over the 
X-intervals where the curves coUapse (0.06 < X < 0.14 at all 
(j>). As a guide to the eye, we also plot the power-law depen- 
dence of s~^° (dashed magenta curve), s"^'^ (dashed violet 
curve) and s~^'^ (dotted cyan curve). Only non-percolating 
clusters are considered in the analysis. The vertical lines 
s* = 50 ± 10 indicates the crossover between two regimes 
at (A = 0.54 . 



exponent t ^ 1.7 for the entire size range (even for very 
small cluster sizes). We notice that this value for r is 
similar to that of the exponent found for the size distri- 
bution of mobile regions in glass- forming systems [60|. 
The discrepancy with the exponent expected from ran- 
dom percolation theory (t '^ 2.2) may stem from the 
fact that new crystalline regions preferably appear in the 
surroundings of existing ones [33] , which causes a partial 
loss of randomness of the growing aggregate. 

We now analyse further the behaviour of /(s) for 
(j) = 0.54 by using CNT, where a clear crossover between 
two regimes is present. When s is small, sub-critical 
clusters form and re-dissolve whereas when s is large 
post-critical clusters irreversibly grow. Only the large-s 
regime seems to follow convincingly a power law with ex- 
ponent T ~ 1, whereas the small-s regime starts deviating 
from random aggregation (r ^ 2.2) when s approaches 
s*. According to CNT [53], the free-energy barrier asso- 
ciated to the appearance of size-s clusters is written as 
/3AG(s) = - In {P{s)), where P{s) = n{s)/N is the prob- 
ability to have size-s clusters in the system. Therefore, 
we can use the cluster size distribution n{s)/N (or equiv- 
alently /(s) x Nc/N) to estimate the barrier, as shown 
in Figure [7] for different packing fractions (averaged over 
the range 0.06 <X < 0.14). 

At all packing fractions higher than 0.55 the curves in 
Figure [7] collapse, showing that crystallites starts grow- 
ing in the same way whenever (f> > 0.55: multiple crys- 
talline patches appear and crystallization proceeds by a 
sequence of stochastic micronucleation events, correlated 
in space by emergent dynamic heterogeneity [33|. As 
we have already mentioned, the case when (j) — 0.55 is 
intermediate between 6 > 0.55 and 6 = 0.54. When 
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FIG. 7: -ln(n(s)/7V) versus s for 0.54 < (j) < 0.61. n{s)/N = 
/(s) X Nc/N (where /(s) data are from Figure[S|. The vertical 
lines s* — 50± 10 indicates the crossover between two regimes 
at (A = 0.54 . 



(j) — 0.55 we cannot distinguish between two regimes and 
we interpret this as a consequence of the practical ab- 
sence of a barrier-crossing nucleation event. Whereas 
when (j) = 0.54 we observe two regimes: one characterised 
by pre-critical clusters (s < s*) and one by post-critical 
clusters at large s. Therefore, in the small s regime, 
— \n{n{s)/N) can be interpreted as the free energy of 
cluster formation whose top can be read from the con- 
tinuous vertical lines in Figures 6 and 7, providing an 
estimate for AG ^ 16. This value is consistent with the 
value of the top of the free energy barrier calculated in 
Ref. [i3, \lM fo'" ^ lower packing fraction (cj) ^ 0.535): 
/3AG* ~ 20. This consistency check makes us gain con- 
fidence of our cluster size distribuition analysis. 



Within Classical Nucleation Theory, once we know 
the top of the free energy barrier /3AG* and the size 
of the critical nucleus s*, it is possible to estimate the 
difference in chemical potential between the fluid and 
the solid: ^A^cnt = ^'^f°' ■ The value of A^icnt 
is order-parameter dependent given that, contrary to 
(3 AG*, s* depends on the choice of the order param- 
eter used to identify solid-like particles IJ]. Only for 
particular choices of the order parameter does A^cnt 
coincide with the true value of A/i (obtained, for in- 
stance, from thermodynamic integration). In our case, 
if we use the CNT expression for (3AficNT taking the 
value of the top of the free-energy barrier and the critical 
nucleus size from Figure [7l^/3AG* — 16 and s* — 50), 
we obtain j3A^cnt — 0.64, which is in good agreement 
with the value computed from thermodynamic integra- 
tion /3A/Z = 0.63 in Ref. lj|. Thus, our choice of the 



order parameter to identify crystalline particles seems to 
be reasonable. 
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D. Effect on crystallization of the cutoff needed to 
define of a solid particle 

We now study the effect of the choice of ^c (the param- 
eter used to define solid-hke particles in our definition 
of crystalhne clusters), on the properties of the clusters, 
such as the radius of gyration. We focus on the two ex- 
treme cases studied, = 0.54 and = 0.61. To identify 
the neighbours of each particle, we use the criterion of 
Ref. [57| and determine the number of connections with 
dg > 0.7. At this point, we use various definitions of the 
number of connections needed for a particle to be identi- 
fied as solid-like: not only ^c = 6 (what we have used so 
far to identify a solid- like particle), but also 8 or 10. To 
qualitatively understand the effect of the choice of ^c on 
the structure of the growing clusters, we represent a slab 
of the system taken at the same position in the sample 
when = 0.54 (Figure [8] a and b) and when = 0.61 
(Figure [5] c and d). In panels a) and c) solid particles are 
evaluated using ^c = 6 and correspond to a system with 
X ~ 0.20 whereas in panels b) and d) they are calculated 
using ^c = 10 and correspond to a system with X ~ 0.10. 




FIG. 8: Slabs (about 3-4 particles diameters thick) where 
solid particles are identified with: Top: = 0.54, a) ^c = 6 
{X ~ 0.20) and b) ^c = 10 (X ~ 0.10); Bottom: = 0.61, c) 
Cc = 6 (X ~ 0.20) and d) ^c = 10 (X ~ 0.10). 

From Figure [51 it is clear that the overall structure of 
the clusters is not much affected by the choice of ^c ■ Fd- 
ion et al. [ij] have recently evaluated the effect of the 
choice of ^c on the calculation of the top of the nucle- 
ation free-energy barrier AG* at = 0.535 (see Figure 
2 in Ref. [IJ]), showing that AG* does not depend on fc 
(named in the same way in their work). This contrasts 
with the corresponding value of the size of the critical 
cluster s*: for ^c = 6 they find s* ^ 100, whereas for 



fc = 10 they obtain s* ^ 25. The authors conclude that 
the main difference among the order parameters (each de- 
fined with a given value of ^c, 5 < ^c < 10) is the ability 
to distinguish between fluid-like and solid-like particles 
near the fluid-solid interface. 

To quantify our observations, we represent the radius 
of gyration as a function of the cluster size at packing 
fraction = 0.54 and = 0.61 (Figure ^ for different 
values of S,c- 
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FIG. 9: Rg versus s at = 0.54 (red) and = 0.61 (black) 
averaged over time and over 10 runs. In both cases, filled 
circles, striped squares and empty diamonds correspond to 
clusters where solid-particles are defined using ^c = &, £,c = 8, 
and ^c = 10, respectively. The spreading of the data for large 
values of s at = 0.61 is due to poor statistics for very large 
clusters. 

Results for — 0.54 and = 0.61 clearly provide evi- 
dence that Rg does not depend strongly on the different 
values of ^c used in the definition of solid-like particles. 
The slabs of Figure |8] corroborate the result presented 
in Figure [HI and we conclude that the structure of the 
growing crystalline clusters is not affected by the choice 

ofCc. 

We now re-plot Figure [8] assigning different colours to 
each particle depending whether a particle has been la- 
belled as solid-like with one of the above-mentioned cri- 
teria: in red, we colour particles that are solid-like ac- 
cording to the most stringent criterion (^c = 10), in or- 
ange particles identified as solid-like using ^c = 8 (but 
not with ^c — 10) and in yellow particles identified as 
sohd-like using ^c = 6 (but not with ^^ = 8). 

When = 0.54 (Figure fTOTleft)) we observe that the 
criterion with ^c = 10 allows one to obtain particles 
in the inner part of the clusters, surrounded by parti- 
cles defined as solid-like according to the looser criterion 
^c = 8 and then by the ones defined with ^c = 6. During 
the time-evolution of the clusters, we always observe the 
colour-code ranging from red (black) in the inside to yel- 
low (light grey) in the outside of the clusters. As we have 
already discussed, small clusters for — 0.54 as well as 
branched ones for = 0.61 are far from being compact: 





FIG. 10: Left: <j) = 0.54. Re-plot of Figure EJi), where we 
colour in red particles identified as solid-like with ^c ~ 10 
(corresponding to the ones represented in Figure [SJj), in or- 
ange particles identified as solid-like using £,c ~ 8 (but not 
with ^c = 10) and in yellow particles identified as solid-like 
using ^c = 6 (but not with ^c = 8). Right: (f) = 0.61. Re-plot 
of Figure [Ht) with the same colour-code. Particles identified 
as solid-like with ^c ~ 10 correspond to the ones represented 
in Figure [Hli). The colours become black, dark grey and light 
grey, respectively, in black-and-white. 



in these cases, most particles will live at the interfaces, 
so that red (black) particles are less abundant. 

In recent works [20-22], it has been suggested that in 
hard sphere systems at low packing fractions (0 ^ 0.54) 
nucleation happens via a two-step scenario: first the for- 
mation of a locally denser regions with high bond orien- 
tational order parameter, and next the restructuring of 
these regions at constant density. Analogously, either at 
(j) = 0.54 and at <j) — 0.61, we first observe the formation 
of loosely packed clusters (in yellow) , that become more 
compact when growing (in orange) and within which we 
identify the growth of highly ordered regions (in red). 
However, from our analysis one could argue that nucle- 
ation happens via a "multi-step" mechanism, and that 
the number of steps depends on the number of cutoffs ^c 
that are monitored, i.e. three in this case. In our previ- 
ous work [331, we have demonstrated how crystallization 
from a glass (at <j) = 0.61) happens via a gradual (many- 
steps) mechanism, but we did not discuss the process at 
packing fractions below the glass transition {(j>g = 0.58). 
As we have just shown, the mechanism of gradual growth 
seems to hold also in super-saturated systems just above 
freezing. 



E. Discussion and conclusions 

In this paper we have analysed the properties of clus- 
ters of solid particles during crystal growth in a system 
of monodisperse hard spheres. We have focused on sev- 
eral packing fractions beyond the fiuid-crystal transition 
up to state points deeply in the glassy regime. Given 
the marked differences in the nucleation processes taking 
place at different (f>, i.e. from CNT regime up to crystal- 
lization from a glass [13, |3j| , we expected to observe clear 
differences also in the properties of the growing clusters. 



On the basis of the results presented here, we can clearly 
identify two regimes where clusters are characterised by 
distinct statistical properties: one pertains to < 0.54 
and one to > 0.56. The value cj) — 0.55 marks the 
threshold between the two behaviours, showing interme- 
diate properties. 

In the \ow-(j) regime {(j) = 0.54), CNT holds and the 
existence of a critical size s* leads to a cluster size distri- 
bution which shows two distinct behaviours respectively 
below and above s* (at long enough times). The struc- 
ture of the clusters is loose for very small sizes (s < 50) 
and crosses to spherical and compact (d/ = 3) for larger 
sizes (this can also be visualised in the snapshots shown 
in Figure[T|). When studying the cluster size distribution, 
we find analogies between crystalline clusters growing 
and standard clusters of attractive particles approach- 
ing percolation: the cluster size distribution seems to 
approach a power-law behaviour also in our case, albeit 
with markedly different exponents. At </) = 0.54 and for 
s < s* , the cluster distribution it is not too different 
from a random-percolation power-law, i.e. /(s) ~ s~^'^. 
However, for s > s* , f{s) crosses to a clearer power-law 
behaviour, with a much lower effective power-law expo- 
nent, i.e. /(s) ^ s^^°, due to the fact that larger clusters 
grow in a compact way. The use of CNT has allowed us 
to provide estimates of the free-energy barrier at ^ = 0.54 
and of the difference in chemical potential between solid 
and liquid phase that, being in agreement with thermo- 
dynamic integration calculations for the same system at 
the same conditions [IJ, ll3| , corroborates our choice of 
the order parameter used to identify solid-like particles. 

In the large-(/) regime {(f) > 0.56), the free-energy bar- 
rier to nucleation becomes negligible, so that no criti- 
cal cluster size exists: small size clusters have a loose 
structure and do not reach a compact form when grow- 
ing {df ^ 2.5). When studying the cluster size distri- 
bution, we observe a power-law behaviour for all sizes, 
but the exponent (r ~ 1.7) is smaller than the value pre- 
dicted by random percolation (r ~ 2.2). This could be 
ascribed to the fact that clusters do not grow completely 
at random, but preferably particles become solid-like in 
the vicinity of solid regions, which act as seeds for further 
growth |33j . This mechanism should become more and 
more important with increasing (jj, due to the reduced 
diffusion of the particles. However, the surprising result 
is that the same behaviour is observed for any packing 
fractions above ((> = 0.56, covering a range where the dy- 
namics slows down by several orders of magnitude [31| . 
Hence, we conclude that the static properties of these 
clusters are not much affected by dynamics. 

We also notice a striking similarity between percolation 
events and nucleation. When (p < 0.56, clusters are more 
compact and take a long time to percolate: similarly, the 
nucleation time (defined as the time when X ~ 0.20 in 
Ref. [l3[) is large at ^ = 0.54, given that crystallization is 
an activate process, and decreases when (j) = 0.55, due to 
the lowering of the free-energy barrier. When (p > 0.56, 
clusters percolate throughout the system and the nucle- 
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ation time increases with (j) due to the slowing down of 
particles dynamics. Of course, the resulting macroscopic 
crystals at all packing fractions are quite different in the 
two nucleation regimes, as evident from the snapshots 
presented in Figure [TJ On the one hand, homogeneous 
and compact structures are found in the CNT region, 
where crystallization once activated proceeds to the for- 
mation of a full crystal. On the other hand, for larger 
(j), the macroscopic crystals that are formed are more 
heterogeneous, since they result from the branching of 
many sub-units. In this case, it is most likely that the fi- 
nal structure will have a poly-crystalline character (with 
different units possibly bearing a different orientation). 

We stress that, for both regimes (low and high 0), the 
cluster properties are very different from those observed 
in a standard aggregating system, driven by attractive 
interactions. The main difference can be observed for 
small sizes, where attraction drives the growth of com- 
pact small clusters (df = 3) [4l|. Then, with increas- 
ing size, the structure of the clusters remains spherical 
upon growth if particles are interacting via a pure attrac- 
tion, while a fractal structure (decreasing df) is observed 
when additional interactions, like for example long-range 
electrostatic repulsion, are also at work [4l|, \^ . In the 
present case of crystalline clusters, the loose structure for 
the small clusters can be understood by the absence of a 
real attraction between the particles. Hence larger clus- 
ters become compact in the classical nucleation regime 
(at small supersaturations) , while they turn into random 
fractal objects when nucleation happens without activa- 
tion (at large supersaturations (j) ^ 0.56), so that adja- 
cent solid particles merge and branch into a percolating 
cluster. 

In conclusion, we have provided new insights into the 
growth of clusters of solid particles during the nucle- 
ation process of monodisperse hard spheres. Our results 
shed more light on the existence of different nucleation 
regimes already discussed in Ref. [Ij, l3ll - [33| , highlight- 
ing the influence of these onto the connective proper- 
ties of the system. In particular we identify that while 
compact crystals are formed at small supersaturations, a 
crossover to the formation of fractal ramified clusters is 
encountered upon increasing (f). Surprisingly, the static 
properties of the growing clusters in the high-0 re gim e 
{(j> > 0.56) are ^-independent. In our previous work [33j, 
we have demonstrated how crystallization from a glass 
(at (f) = 0.61) happens via a gradual (many-steps) mech- 
anism, but we did not discuss the process at packing 
fractions below the glass transition {(pg — 0.58). In this 
paper we show how the mechanism of gradual growth 
seems to hold also in super-saturated systems just above 
freezing. 

The two growth regimes (at < 0.56 and (/> > 0.56) 
may explain the different morphology of the crystals ob- 
served at low and high cj) in the experiments of Pusey 
and van Megen [2|. However, to allow for a more detailed 
comparison with those experimental results, a similar in- 
vestigation for slightly polydisperse hard spheres, which 



should also favour fractionation and hence a stronger ten- 
dency to the formation of poly-crystals, will be the sub- 
ject of a future investigation. 
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Appendix A: Local bond-order parameter 



We use the 
der parameter 
tween particle 
each of its neighbour j, d(,[i,j) 



rotationally invariant local bond or- 
dg defined as the scalar product be- 
i's gg complex vector and the one of 

El=-6 96,m(«) • 



le miJ) (where q^ ^ is the complex conjugate), then 
averaged over all particle i's neighbours, Nb{i) pU - 
|63[. Each component qe^mii) depends on the rela- 
tive orientation of particle i with respect to its Ni,{i) 
neighbouring particles, and is defined as qe.mii) ~ 

)] I teL=-6 k6,™(*)P (with 



Wb(i) 



are the spherical harnion- 



X] j=l V'e.m (^i J- J Y'J J ; 

TO = [-6,6]), where -06,m( 
ics of order 6. 

In order to identify particle «'s neighbours we use the 
criterion of Ref. [531, that is capable of avoiding a den- 
sity dependence cutoff distance r^. We have checked that 
such criterion is consistent with a cutoff distance crite- 
rion, where r\, is tuned for any packing fraction. We then 
consider particles i and j as having a " solid connection" 
when their d%(i,2) exceeds the value of 0.7: particle i 
is then labelled as solid-like if it has at least ^c=6 solid 
connections. 

dQ{i,j) is a normalised quantity correlating the local 
environments of neighbouring particles, it is a real num- 
ber and is defined in the range —1 < de{i,j) < 1: it 
decreases when thermal vibrations are present but, on 
average, it is close to one if particles have a solid-like en- 
vironment, and around zero if particles have a fluid- like 
environment. 



Appendix B: Values of Ast used to compute /(s) 



In section IIII CI we compute the fraction of clusters of 
size s, /(s): each value of /(s) is computed within a given 
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interval of width Asi = {si,max ~ Si^min)- In Table |T] we 
explicitly indicate all Asi used to calculate /(s). In our 
calculations we have also checked that different choices of 
interval widths yield the same cluster size distributions, 
within statistical error. 



TABLE I: Bonudaries of the cluster size intervals, 



Ji.max , -31,771 



Si,m.in], centred at s. 



[Si^max: Si,min 


s 


\Si,7naxi Si^rnin 


s 


\Si,maxt Si^rnin\ 


s 


1 


1 


[36,45] 


41 


[451,600] 


525 


2 


2 


[46,55] 


53 


[601,1000] 


800 


3 


3 


[56,75] 


65 


[1001,1300] 


1150 


[4,7] 


5 


[76,95] 


85 


[1301,2000] 


1650 


[8,11] 


9 


[96,120] 


105 


[2001,3500] 


2750 


[12,15] 


13 


[121,140] 


130 


[3501,5000] 


4650 


[16,21] 


17 


[141,200] 


170 


[5001,7000] 


6000 


[22,25] 


23 


[201,250] 


225 


[7001,10000] 


8500 


[26,27] 


26 


[251,300] 


275 


[10001,15000] 


12500 


[28,35] 


31 


[301,450] 


375 
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